home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
DDJMAG
/
DDJ9110.ZIP
/
EIGENCAL.ASC
< prev
next >
Wrap
Text File
|
1991-09-10
|
9KB
|
310 lines
_MIXED-LANGUAGE WINDOWS PROGRAMMING_
by John Norwood
[LISTING ONE]
subroutine tred2(a,n,np,d,e)
c
c Householder reduction of a real symmetric, nxn matrix a, stored in an
c npxnp physical array. On output, a is replaced by the orthogonal matrix
c q effecting the transformation.
c d returns the diagonal elements of the tridiagonal matrix, and e the off-
c diagonal elements, with e(1)=0.
implicit none
integer*4 i,j,k,l ! Loop indices
integer*4 n ! Actual array size
integer*4 np ! Physical array size of incoming array
real*4 a(np,np) ! Matrix, np is used so dimensions match
real*4 d(np) ! Vector that will have diagonal elements
real*4 e(np) ! Vector that will have off-diagonal elements
real*4 h ! Vector norm used in forming projection
real*4 scale ! Scale factor
real*4 f ! Temporary variable
real*4 g ! Temporary variable
real*4 hh ! Another piece of the projection
if (n.gt.1) then
do 18 i=n,2,-1
l=i-1
h=0.
scale=0.
if(l.gt.1) then
do 11, k=1,l
scale=scale+abs(a(i,k))
11 continue
if(scale.eq.0.) then ! Skip transformation
e(i)=a(i,l)
else
do 12, k=1,l
a(i,k)=a(i,k)/scale ! Use scaled a's in transformation
h=h+a(i,k)**2 ! for eigenvectors
12 continue
f=a(i,l)
g=-sign(sqrt(h),f)
e(i)=scale*g
h=h-(f*g)
a(i,l)=f-g
f=0.
do 15,j=1,l
a(j,i)=a(i,j)/h
g=0.
do 13, k=1,j
g=g+a(j,k)*a(i,k)
13 continue
if(l.gt.j) then
do 14,k=j+1,l
g=g+a(k,j)*a(i,k)
14 continue
endif
e(j)=g/h ! Form element of projection in
f=f+e(j)*a(i,j) ! temporarily unused element of e
15 continue
hh=f/(h+h)
do 17, j=1,l
f=a(i,j)
g=e(j)-hh*f
e(j)=g
do 16, k=1,j ! Loop to reduce matrix a
a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
16 continue
17 continue
endif
else
e(i)=a(i,l)
endif
d(i)=h
18 continue
c This starts the eigenvector specific part of the code.
endif
d(1)=0.
e(1)=0.
do 23, i=1,n ! Begin accumulation of transformation matrices
l=i-1
if(d(i).ne.0.) then
do 21, j=1,l
g=0.
do 19, k=1,l ! Use information stored in a to form
g=g+a(i,k)*a(k,j) ! projection times orthogonal matrix
19 continue
do 20, k=1,l
a(k,j)=a(k,j)-g*a(k,i)
20 continue
21 continue
endif
c This ends the eigenvector specific part of the code.
d(i)=a(i,i)
c This starts the eigenvector specific part of the code.
a(i,i)=1. ! Reset row and column of a to identity matrix
if(l.ge.1)then ! for the next iteration of transformation loop
do 22,j=1,l
a(i,j)=0.
a(j,i)=0.
22 continue
endif
c This ends the eigenvector specific part of the code.
23 continue
return
end
subroutine tqli(d,e,n,np,z,iterflag)
c This subroutine performs a QL algorithm with implicit shifts, to determine
c the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix,
c or of a real, symmetric matrix previously reduced by tred2 above.
c d is a vector of length np. On input, its first n elements are the
c diagonal elements of the tridiagonal matrix. On output, it returns the
c eigenvalues. The vector e inputs the subdiagonal elements of the
c tridiagonal matrix, with e(1) arbitrary. On output, e is destroyed.
c If eigenvectors are desired, the matrix z (nxn stored in an npxnp array)
c is input as the identity matrix or the matrix that is returned from tred2.
c On output, the kth column of z contains the normalized eigenvector
c corresponding to the eigenvalue in d(k).
implicit none
integer*4 i,k,l ! Loop indices
integer*4 iter ! Iteration counter
integer*4 n ! Logical size of array z
integer*4 m ! Submatrix size
integer*4 np ! Physical size of array z
integer*4 iterflag ! Flag to return error code to main routine
real*4 d(np) ! Diagonal elements is, eigenvalues out
real*4 e(np) ! Off diagonal elements in, nothing out
real*4 z(np,np) ! Matrix from tred2 in, eigenvectors out
real*4 dd ! Holds small subdiagonal element
real*4 g ! Holds Givens rotation
real*4 s ! "Sin" component of Givens rotation
real*4 c ! "Cos" component of Givens rotation
real*4 p ! Element of projection matrix
real*4 f ! Holds result of "sin" applied to element of e
real*4 b ! Holds result of "cos" applied to element of e
real*4 r ! Temporary piece of c or s
if(n.gt.1) then
do 11, i=2,n ! Renumber elements of e
e(i-1)=e(i)
11 continue
e(n)=0.
do 15,l=1,n
iter=0
1 do 12,m=l,n-1 ! Search for small subdiagonal element
dd=abs(d(m))+abs(d(m+1))
if (abs(e(m))+dd.eq.dd) go to 2
12 continue
m=n
2 if(m.ne.l) then
if(iter.eq.30) then
iterflag = -1
return
endif
iter=iter+1
g=(d(l+1)-d(l))/(2.*e(l)) ! Calculate shift
r=sqrt(g**2+1.)
g=d(m)-d(l)+e(l)/(g+sign(r,g))
s=1.
c=1.
p=0.
do 14, i=m-1,l,-1 ! Plane rotation followed by a Givens
f=s*e(i) ! rotations to maintain tridiagonal
b=c*e(i) ! form
if(abs(f).ge.abs(g))then
c=g/f
r=sqrt(c**2+1.)
e(i+1)=f*r
s=1./r
c=c*s
else
s=f/g
r=sqrt(s**2+1.)
e(i+1)=g*r
c=1./r
s=s*c
endif
g=d(i+1)-p
r=(d(i)-g)*s+2.*c*b
p=s*r
d(i+1)=g+p
g=c*r-b
c Start of code specific to forming eigenvectors.
do 13, k=1,n
f=z(k,i+1)
z(k,i+1)=s*z(k,i)+c*f
z(k,i)=c*z(k,i)-s*f
13 continue
c End of code specific to forming eigenvectors.
14 continue
d(l)=d(l)-p
e(l)=g
e(m)=0.
go to 1
endif
15 continue
endif
return
end
Examplσ 1:
(a)
SUBROUTINE STRINGER(INSTRING)
CHARACTER*40 INSTRING
INSTRING = 'This is from Fortran'
END
(b)
è
Sub Form_Click ()
Dim temp As String * 40
Call STRINGER(temp)
Debug.Print temp
End Sub
(c)
Declarσ SuΓ STRINGE╥ LiΓ "d:\vb\test\string.dlló (ByVa∞ Mystrinτ A≤ String⌐
Examplσ 2:
(a)
SUBROUTINE ARRAYTEST(ARR)
INTEGER*4 ARR(20)
ARR = 5
END
(b)
Sub Form_Click ()
Static testarray(1 To 20) As Long
Call ARRAYTEST(testarray(1))
For i% = 1 To 20
Debug.Print testarray(i%)
Next i%
End Sub
(c)
Declarσ SuΓ ARRAYTES╘ LiΓ "d:\vb\test\array.dlló (Myarra∙ A≤ Long)
Examplσ 3
(a)
SUBROUTINE ARRAYSTRINGS(ARR)
CHARACTER*24 ARR(5)
ARR = 'This is a string also'
END
(b)
Sub Form_Click ()
Static testarray(1 To 5) As StringArray
Call ARRAYTEST(testarray(1))
For i% = 1 To 5
Debug.Print testarray(i%).strings
Next i%
End Sub
(c)
Type StringArray
strings As String * 24
End Type
Declare Sub ARRAYSTRINGS Lib "d:\vb\test\array.dll" (Myarray As StringArray)